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В статье предлагаются новые коллокационные блочные методы решения жестких систем обыкновенных 
дифференциальных уравнений с автоматическим управлением размером шага, которые изначально 
ориентированы на параллельную архитектуру. Основная идея, на которой базируется конструирование 
предлагаемых методов, заключается в одновременном получении приближений точного решения во 
всех расчетных точках блока. 
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порядок аппроксимации. 
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Введение 


Основные усилия при параллельном моделировании динамических объектов 
большой размерности, описываемых системами обыкновенных дифференциальных 
уравнений (СОДУ), в настоящее время концентрируются на распараллеливании 
последовательных алгоритмов и программ [1], [2]. 
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При этом исследуются только информационные взаимодействия смежных час- 
тей алгоритма, и не уделяется должного внимания общей информационной структуре 
задачи, что не позволяет создать масштабируемые параллельные программы [3], [4]. 

Еще одной проблемой при моделировании динамических объектов, особенно 
тех, которые формализуются жесткими системами ОДУ, является необходимость 
изменения шага интегрирования [5], [6]. 

Также актуальными являются вопросы соотнесения погрешностей результатов 
и времени, затрачиваемого на получение решения. 

Практически все известные в настоящее время численные методы с автома- 
тическим выбором шага интегрирования основаны на вычислении главного члена 
локальной ошибки и последующем выборе такого размера для очередного шага, 
который является максимальным для заданного предела локальной точности [5-8]. 

Но при этом возникает необходимость повторных вычислений в каждой точке, 
что приводит к значительному увеличению вычислительных затрат. 

Кроме того, жестко регламентируются пропорции сокращения шага [9], [10]. 

В статье предлагаются новые блочные методы решения систем обыкновенных 
дифференциальных уравнений, которые изначально ориентированы на параллельную 
архитектуру [11]. 

Основная идея, на которой базируется конструирование блочных методов для 
решения СОДУ на параллельных компьютерах, заключается в одновременном получе- 
нии приближений точного решения в точках, образующих блок. 

Поскольку точки внутри блока расположены регулярно, есть возможность опреде- 
ления и сопоставления локальных погрешностей во всех точках блока, в отличие от 
стадийных методов, где стадии, как правило, не совпадают, и сравнение ведется только 
по одной конечной расчетной точке 1,„.у. 


Также непривлекательным является использование в компьютерах с распределен- 
ной памятью мелкозернистого параллелизма стадийных методов. 

Алгоритмы управления шагом, предлагаемые в статье, базируются на исполь- 
зовании коллокационных блочных одношаговых и многошаговых методов. Парал- 
лельный счет осуществляется в пределах одного цикла для всех точек блоков с 
размерностями 5 и 5+ Г. 

Две нити вычислений проходят независимо, а необходимость в обменах возникает 
только после получения конечных результатов для блоков расчетных точек. 

Цель данной статьи состоит в создании коллокационных блочных расчетных 
схем, позволяющих обеспечивать автоматическое управление шагом интегрирования 
при моделировании динамических объектов в параллельных вычислительных системах. 


Генерация коэффициентов расчетных 
схем коллокационных блочных методов 
Для параллельной реализации численного решения задачи Коши 
= У(в( же), Жю)=м (О 


можно обобщить подходы, предложенные в [12], [13], и основанные на одношаговом 
и многошаговом блочных методах. При этом не обязательным является требование, 
связанное с совпадением размерностей опорного и расчетного блоков. Будем счи- 
тать, что схема одношагового коллокационного блочного метода содержит 5 узлов в 
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рассчитываемом блоке (рис. 1) и одну опорную точку, многошагового -— 5 узлов в 
рассчитываемом блоке при использовании вычисленных значений приближенного 
решения в т предшествующих блоку точках (рис. 2). 


[-15=[,0 [1 [2 В5=Ё-10 [+11 


у 


блоки 


Рисунок 1 — Шаблон разностной схемы одношагового 
5-точечного блочного метода 


ВЕ=Ь 1,0 


Г и,-(т-1) 


точки опорного (опорных) блоков точки расчетного блока п 


Рисунок 2 — Шаблон разностной схемы т-шагового 
5-точечного коллокационного блочного метода 


Уравнения многошаговых разностных методов для блока из 5 точек при 
использовании вычисленных значений приближенного решения В т 


предшествующих блоку узлах, с учетом введенных выше обозначений, можно 
записать в виде: 


> аш,,= ей. == [,2,....5, п= 1,2... (2) 


1=1-т 1=1-т 


- 
т 


Для одношаговых методов, представляющих частный случай (2), если т=[, 
разностные уравнения имеют вид 
] 5 5 : 
Е. От: ан ле чае Е а п 


Формулы (3) определяют одношаговый 5-точечный разностный метод. 
Обозначим матрицы коэффициентов системы уравнений (2) через 


А={а,,р В={Ь,,1=1,2,...5, у=-(т-1)-(т-2),...,5. 
Разобъем каждую из матриц на две части 
А ={а;, В, 1=1,2,...5, ЛД =-(т-1),-(т-2),...,0. 
А» ={а, ;}, В? ={Ь;; и У 
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Введем соответствующие вектора 
Ив={и,,/) Л=-(т-1)-(т-2),...,0, 
У, 1 Инн Я 
Ри ={ ЛА 1=-(т-1),-(т-2),...0 


О 
В векторной форме уравнение (2) будет иметь вид 


и о о ми ВО. сы (4) 


Чтобы система однородных разностных уравнений, соответствующая (4), была 
линейно независима, потребуем, чтобы матрица А› была невырожденной, и разрешим 


ее относительно У», } 


=О0, +т(РЕ, +СЕ,), (5) 


У 


где О=-АА,Р=А В, С=АВ.. 

Задав стартовые значения, полученные, например, с помощью явных стадийных 
схем, необходимо на каждом последующем этапе решить нелинейное уравнение (5), 
определив последовательно вектора У), \У»,..... 

Представление блочных разностных многошаговых многоточечных уравнений 
в виде (5) будем называть канонической формой их записи. Получим коэффициенты 
уравнений, позволяющие представлять методы (2) в канонической форме (5). Каждое 
1-е разностное уравнение 


а кт авы №: 1=12,... (6) 


1=1-т 1=1-т 


содержит 271+5 неизвестных коэффициентов 4;,,А„;, = (т-1)—т-2),....0 
8„» /=1,2,....5. Для их определения следует использовать 27+5 уравнений условий 


аппроксимации. Выражения для невязок на решении х(1) исходного дифференциаль- 
ного уравнения имеют вид 


а т а а. м изв, (7) 


]=-т ДЕ-т 


где Х„; = (1, на Л ), Хх п,0 } ХЕХЕ, =а Л )= (Е, че д ), 


' ' 


=х 


=х 


п-т 


Хх 


п-1,т п,0 * 


Для каждого 1-го уравнения потребуем его аппроксимации в точке #, 


Зе Е о+ п) Хи ЕХ(Ьо+ Л = Уч, ли ). 
Раскладывая х„; = (0+ М) их, =Х(1 0+ Л) в ряды Тейлора в окрест- 
ности точки 1,0), подставляя эти разложения в выражение (7), группируя члены с 


одинаковыми степенями по т, а затем приравнивая нулю коэффициенты при т’, 
получим систему уравнений 
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0 
а о: 


1=1-т 


0 5 
>04: +4,;)+ вы =, 
= 


1=1-т 


0 .[ $ .[ 
] 1-1 аа 
Я а НЕ = Е (8) 
›] [ >] >] [ 
]Д=Е-т 9=1 

Поскольку каждое 1-е разностное уравнение (6) содержит 2т+5 неизвестных 
коэффициентов, то максимальный порядок аппроксимации рассматриваемого т-ша- 
гового $-точечного разностного метода равен р=2т+ 5-1. 

Каноническая форма одношаговых многоточечных уравнений будет иметь 
соответственно вид 


5 
пин о (оо + Утв, (9) 

1=1 
Коэффициенты разностных уравнений (6) и (9) можно определить, решая 
систему (8) для общих многошаговых блочных методов, учитывая их частный вид, 
или интегро-интерполяционным методом. Для этого строится интерполяционный 


многочлен Г„.,_/(1) с узлами интерполяции 1, ,_„ и соответствующими им значе- 


ниями сеточной функции ии, 1= -(т - 1),-(т-2),...,0,1,...5. Находятся производные 


полученного интерполяционного многочлена в узлах сетки [,;,1= [,2,...,5. Приравняв 


п? 
их соответствующим значениям правой части дифференциального уравнения, получим 
разностные уравнения для блока. Используя изложенный выше подход, можно 
сформировать каноническую систему уравнений для любых значений параметров т 
и 5. В частности, для параметров т=3, 5=2, р=7 и шаблона, приведенного на рис. 3, 
система будет выглядеть следующим образом 


аи ра га м ой Я То 
Ул, ых 1360 170 170 У 170 1360 № 
Ул, 36}. Р-Я 806 1,0 44 Тел Г 247,2 
85 85 85 85 85 
ЭЙ. 5 ИИ ре ЛЭСИ 
, рт Сер ЗАВ з 
4: 136 17 136 
54и,_› 1281 т 199и 0 
17 17 17 


блок и-1 блок и 


Рисунок 3 — Шаблон разностной схемы коллокационного 3З-шагового 
2-точечного блочного метода 
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Управление шагом при интегрировании 
блочными методами 


Рассмотрим решение задачи (1) одношаговым коллокационным блочным 
методом (9) с числом расчетных точек 5. Выделим две системы процессорных узлов, 
на которых запустим параллельные процессы. Первая система узлов будет 
обеспечивать параллельную реализацию одношагового коллокационного блочного 5-то- 
чечного метода, что потребует 5 процессорных элементов для получения значений 


Ии+ 11, [= [,2,...,5. На второй системе узлов будет осуществляться реализация одно- 


п- 


шагового коллокационного блочного 5+/-точечного метода с получением значений 
-1;, = 1,2,...5+1. Поскольку методы являются неявными, каждый временной 


ил 
шаг подразумевает проведение некоторого количества итераций, обеспечивающих 
требуемую локальную точность. 

Реализация обоих методов будет осуществляться автономно, что обеспечит 
крупнозернистость вычислений. 

Размерность вычислительного поля может быть представлена двумя вариантами. 
В первом случае это могут быть кольцевые топологии, количество процессорных 
элементов в которых совпадает с размерностью блоков 5 и 5+/ соответственно, во 
втором случае вычислительные поля представляют собой решетки процессорных 
элементов (рис. 4). 

Первая решетка с 5 (количество точек в блоке) столбцами и М (количество 
уравнений в системе) строками, у второй решетки размерность процессорного поля — 
(5+1)хМ. 

Если используется кольцевое вычислительное поле, размерность которого 
совпадает с размерностью блока, каждый процессорный элемент закрепляется за 
рассчитываемой точкой блока для всех уравнений системы. 

Для осуществления итераций по (9) в каждом 1-м процессорном элементе, 
обеспечивающем вычисления 5-точечным методом, должны быть размещены 
соответствующие коэффициенты а; ;,Б;, [= 1,2,...,5, ] = 1,2,...5, а также значения 
элементов вектора правых частей системы в последней точке предшествующего 
блока, которые будут считаться начальными для следующего. 

Для группы процессорных элементов, которые обеспечивают решение 5+ /-то- 
чечным методом, соответствующие коэффициенты а; ‚,Б;, {= 1,2,...,5+1, } = 1,2,....5+1. 


В случае, если используется решетка процессорных элементов, позиция каж- 
дого элемента в строке определяет порядковый номер уравнения в системе, для 
которого ведется расчет, а номер элемента в столбце определяет рассчитываемую 
точку блока. 

Поскольку вычисления проводятся для блоков точек, которые расположены 
регулярно, есть возможность сопоставления решений, полученных методами с раз- 


ными порядками точности в 5-совпадающих расчетных точках 1„.1,1„-2,..., .;. ЕСЛИ 
норма вектора расхождений не превосходит заданную глобальную точность 


<ю[,1 = [,2,....5, за основу берется решение, полученное 


= 1[,2,....5+Г[. 


вычислений | и 


м Уы 


коллокационным блочным 5+ /-точечным методом У „.;, 


Рассчитывается новое значение шага т 


пем 
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Слеи = Ит| Дахтах, тах Гахпит, © 


и по результатам расчета 5+ /-точечным методом формируется новый вектор опор- 


ных точек. 
Е ы Параллельное вычисление 
п+19 э Чи ли- 19+ * М п-+1.0 многошаговыми 
коллокационными блочными 
методами с числом точек $ и $+1 


Ти-1ь 5 Инн и-1° .. °И+1,0 


. 
«- 


Формирование нормы вектора 
расхождений 
в совпадающих точках 


да Формирование решения 


Сокращение по блоку 
шага Ти+1 размерностью $+1 


нет Хин Уп 


ое 


Хот, 5+1 Уп+1, 51 


нет 


Формирование решения 
по блоку размерностью $ 


Хи-+1 1 -— Мп, 1 


...9 


Хп+1,3— Ип+1,5 


Тп+2 » Ии+2,иь- 15+ « «М п+2,0 


Тп+25 $5 Мп+2т-19* + + +2,0 Тп-+2. 5+1 э Шт па- 19° «9 1+2,0 


Формирование нового 
шага и начального 
приближения 


Рисунок 4 — Алгоритм управления шагом интегрирования 
на решетке процессорных элементов 


Если норма вектора расхождений не превосходит заданную локальную 
<=, 1 = [,2,...,5 , за основу берется решение, 


точность вычислений, т.е. 10[ < [м У: 
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полученное коллокационным блочным 5-точечным методом, их, [= 1,2,...,5. Рас- 


считывается новое значение шага т„„,, и новый вектор опорных точек формируется 


по результатам расчета 5-точечным методом. Если полученное решение не обеспе- 
чивает заданную локальную точность, от шага необходимо отказаться, сократив его 
на величину, задаваемую параметром /ахпит. 

Параллельное управление шагом для задачи (1) многошаговым коллокационным 
блочным методом (6) с числом опорных точек т и расчетных 5 не будет иметь 
принципиальных различий с подходами, рассмотренными выше. Однако возникает 
необходимость в формировании начальных данных для расчета очередного блока 
значений. Для выполнения следующего шага нужно обеспечить доступ каждого про- 
цессорного элемента к значениям правых частей в опорных точках блока. 

Кроме того, для осуществления итераций по (6) в каждом 1-м процессорном эле- 
менте, обеспечивающем вычисления 5-точечным т-шаговым методом, должны быть раз- 
мещены соответствующие коэффициенты а; ‚,Б;, {= 1,2,...,5, } = 1,2,...,5,[ = [,2,...т, а 


также значения элементов векторов правых частей системы в т точках предше- 
ствующего блока, которые будуг считаться начальными для следующего. Для группы 
процессорных элементов, которые формируют решение 5+ /-точечным методом, со- 
ответствующие коэффициенты а, ‚В; 1, [= 1,2,...5+1, ] = 1,2,...,5+1,1= 1,2,...т. 


Поскольку максимальный порядок р аппроксимации рассматриваемого т-ша- 
гового 5-точечного разностного метода, который равен р=2т+5-[, не обеспечивает 
абсолютной устойчивости метода [12], при генерации расчетных коэффициентов 
необходимо сократить предельное значение, введя параметр 4<р. Тогда, в зави- 
симости от значения нормы вектора расхождений, будет изменяться расчетная схема 
определения шага интегрирования и новое значение шага *„„, определяется с учетом 


пем 


поправки на порядок 4 аппроксимации. 


Реализация тестовых задач на основе 
алгоритмов управления размером шага 


Для экспериментов на мультиосновных машинах использовался программный 
интерфейс МаГлиК среды параллельных вычислений в Мафетайса \УоШат 
Кезеагсв, обеспечивающий возможность параллельного поиска решений как на 
локальной машине, так и в сетевом масштабе. 

Также параллельная реализация осуществлялась на кластере МеСа5 МТМО- 
архитектуры с распределенной памятью. 

Кластер состоит из 93-х вычислительных узлов, узла управления (Егопё Моде), 
системы коммутации в составе двух гигабитных Еегие{-коммутаторов (З\уйсв НР 
Ргосигуе). 

В качестве тестовой выбрана жесткая задача из [9] 
5(хз 


о ее ба 0 2%, (10) 


№ = ТА, д, Эта 
с начальными условиями 
х(0)=х.(0)=х.(0)=х.(0)=1 
и с известными точными решениями 


&1т( 12 ) 55т( 12 ) 


х(Е)=е ‚Жо(Е)=е ‚ же) = т(Е? )+Ь х.(Е)= с05( 1? ). 
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Поиск решения проводился на интервале [0, 2.5]. Из-за разных масштабов решения, 
особенно в середине интервала интегрирования, потребовалось использовать проце- 
дуру управления шагом, основанную на рассмотренных коллокационных блочных 
методах (рис. 5а). Рассматривались реализации на основе одношаговых коллокационных 
блочных методов с числом расчетных точек 2 и 3, и многошаговых -— с числом опор- 
ных точек 2 и расчетных 2 и 3. Шаги осуществлялись параллельно, продвижение 
вдоль оси времени проходило сразу на размерность блока, обеспечивающего задан- 
ную точность, к середине интервала интегрирования заметно сокращение шага. 

Вторая тестовая задача из [14] 


х’=Ь-х] Хх», х›'= -100х, —х>, (11) 
с начальными условиями 
х(0)=1, х›(0)=0, х;(0)=1, х.(0)=0, 
имеет проблемы при численном интегрировании в начале отрезка (рис. 56) 


ер 
Мер 
0.08 | 1 
0.04 | 1 
0.061 1 0.03 | . н | ы 
т = ый РА ` - * 
0.02 | ТАНИ ИЯ А КВА 
0.011 Л 
фе ии] | 
пер пер 
20 40 60 80 100 120 
6) 


Рисунок 5 — Управление шагом интегрирования для тестовых задач (10) и (11) 


Тестовые жесткие задачи третьего типа [15] возникают при дискретизации 
уравнения диффузии методом прямых. Количество уравнений в системе М может 
быть сколь угодно велико и определяется размерностью расчетной сетки, которая 
вводится для решения исходного уравнения в частных производных. Система (1) 
решалась на интервале [0, 4], правые части системы имеют вид 


-2 1 0 0.00 0 

ый 1 0.00 0 
ям ОО бы 0. |, а 

0000 259. 1 0 

0000 т 59 Ф(#). 


где 
ф(Ё)=@ зат ехр(- 9#)— зж(Т)ехр(- м), 
= а = с05(/2 )/ ыо Сео т) 
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3 = 2(м +1}? 1-с05 


М+ 1 


Точные решения 


1 | 
х(Е)=а т ехр(-= 91 )— 5т ехр(-= иг), 1=1,2,..,М. 
№М+1 
На графиках приведены численные результаты решения (рис. ба) и погрешности 
(рис. 66), полученные для системы с М№ = 8, а также диаграмма параллельного 


управления шагом (рис. 7). 


х1, хо, х3, х4, х5, хб, Х/, х8 


{Етог х1, х2, ХЗ, х4, х5, хб, х7, х8 
Цою 61| 


5 0Ю 71 


5.00 7: 


6) 


Рисунок 6 — Численное решение и погрешности системы (12) 


Мер 
0.20 


0.15 
0.10 
0.05 


озер 


Рисунок 7 — Управление шагом интегрирования для тестовой задачи (12) 


Выводы 


В статье рассмотрены вопросы разработки коллокационных блочных методов 
решения задачи Коши, ориентированных на использование в параллельных вычис- 
лительных системах с распределенной памятью. Основная идея, на которой базировалось 
конструирование блочных методов для решения СОДУ на параллельных компью- 
терах, заключалась в одновременном получении приближений точного решения в 
точках, образующих блок. Параллельный счет осуществлялся в пределах одного цикла 
для всех точек блоков с размерностями 5 и 5+1. Две нити вычислений проходили 
независимо, и необходимость в обменах возникала только после получения конечных 
результатов для блоков расчетных точек. Структура предлагаемых методов позволила 
построить алгоритмы автоматического управления шагом интегрирования, что явля- 
ется особенно актуальным при моделировании объектов, которые описываются жесткими 
или плохо обусловленными системами уравнений. Поскольку точки внутри блока 
расположены регулярно, стало возможным определение и сопоставление локальных 
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погрешностей во всех точках блока, в отличие от стадийных методов, где стадии, как 
правило, не совпадают, и сравнение ведется только по одной конечной расчетной точке. 

На основе предложенных алгоритмов управления размером шага реализованы 
тестовые задачи, численное решение которых обеспечивает заданную точность с 
максимально возможным при этом шагом интегрирования. 

Для экспериментов на мультиосновных машинах использовался программный 
интерфейс МаГлик среды параллельных вычислений в Мафетайса \/оШат Кезеагсв, 
обеспечивающий возможность параллельного поиска решений как на локальной ма- 
шине, так и в сетевом масштабе. 
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оп ше РагаЙе! СоПосанпоп Воск Метоа$5 


Тре рарег 415сиззе Ше 155ие оЁ Фе 4езепт оЁ соПосайоп тешо4$ Гог зо]уше 
Саиспу ргоет, \мсВЬ аге опегие Гог изше ш рагаПе сотрщег$ \уий зе щед 
тетогу. Тре таш 14еа, умсЬ ипдегпез фе деуе]ортепе оЁ фе Моск тефо@$ Гог 
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Фе ехас( зошНоп аё ро Тогтше Фе Ыоск. А рагаПе| соипе \аз регогтед \уифт опе 
субе Гог аП рошб оЁР Моск$ \ий 4итеп$101$ оЁ 5 ап4 5 +1. Туо Но\з оЁ сотршаНоп 
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